Analisi dell’espressione differenziale in diverse condizioni

Nella seguente analisi, l’obiettivo è studiare l’espressione differenziale di specifici geni in diverse condizioni, in particolare nel contesto del cancro al seno. Utilizziamo esemplari di topo, che condividono gran parte del genoma con l’essere umano, per studiare l’espressione differenziale in tre condizioni: silenziamento di BRCA1 (gene che produce una proteina coinvolta nella riparazione del DNA; la sua disfunzione è associata a un aumentato rischio di sviluppare certi tipi di cancro, come il cancro al seno), silenziamento di BRCA2 (gene che produce una proteina fondamentale per la riparazione del DNA e la stabilità genomica. Le proteine codificate da questi geni aiutano a riparare i danni al DNA che possono portare a mutazioni e sviluppo di tumori), e la condizione “Wild Type” (senza silenziamento).

Analisi Silenziamento BRCA2 vs WT senza silenziamento

In questo primo confronto si intende confrontare l’espressione differenziale tra i campioni tumorali che presentano un silenziamento del gene BRCA2 e i campioni tumorali che non presentano nessun tipo di silenziamento

Step 1: Import dei dati

Sfruttiamo la libreria GEOquery per caricare i dati relativi a questo specifico esperimento e successivamente consideriamo solo la matrice dei conteggi grezza di interesse da cui estraiamo i conteggi

library(GEOquery)
gse <- getGEO("GSE137818")[[1]]
#getGEOSuppFiles("GSE137818")
raw_data <- read.csv("GSE137818/GSE137818_Mouse_UNTREATED_bulkRNA_BRCA2KO_vs_WT_rawcounts.csv.gz", stringsAsFactors = FALSE)
head(raw_data[,1:5])
##           X KO_d14_IgG_1 KO_d14_IgG_2 KO_d14_IgG_3 KO_d14_IgG_4
## 1 100009600           23           29           57           28
## 2 100009614            0            0            0            0
## 3 100009664            1            1            0            0
## 4    100017         2330         2586         2683         2705
## 5    100019         3428         3522         2736         3616
## 6 100033459          217          265          166          101
conteggi <- raw_data[,-1]

Step 2: Modifica e processing dei dati

In questa fase andiamo a modificare i nomi dei campioni rispetto al formato orginale per migliorarne l’interpretazione

old_names <- strsplit(names(raw_data), split="_")
new_names <- vector("list", length(old_names))

for (i in 2:length(old_names)) {
  if (old_names[[i]][1] == "KO") {
    old_names[[i]][1] <- "BRCA2"
  }
  if (length(old_names[[i]]) > 1) {
    new_names[[i]] <- paste(old_names[[i]][1], old_names[[i]][length(old_names[[i]])],  sep = "-")
  } else {
    new_names[[i]] <- old_names[[i]][1]
  }
}
new_names <- unlist(new_names)

names(conteggi) <- new_names #abbiamo cambiato i nomi rendendoli un po' più leggibili
head(conteggi)
##   BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2 WT-3 WT-4 WT-5
## 1      23      29      57      28      28   16   24   23   14   10
## 2       0       0       0       0       0    0    0    0    0    0
## 3       1       1       0       0       0    0    0    0    0    0
## 4    2330    2586    2683    2705    3069 1659 1978 1973 1437 1591
## 5    3428    3522    2736    3616    3941 2353 3357 2990 2366 3362
## 6     217     265     166     101     119   36  153   77   59   52

Andiamo a costruire l’oggetto della classe SummarizedExpreriment, ma prima andiamo a costruire i rowData formati dalla coppia ENTREZID-GENE SYMBOL per ogni gene, i colData formati semplicemente dai nomi dei campioni e dalla condizione Silenziato/non silenziato

library(SummarizedExperiment)
library(clusterProfiler)
id_geni <- raw_data[,1]
rownames(conteggi) <- id_geni

organism <- 'org.Mm.eg.db'
gene_symb <- bitr(id_geni, fromType="ENTREZID", toType="SYMBOL",OrgDb = organism)
row_data <- data.frame(ENTREZID=id_geni, ordine=seq_along(id_geni))
row_data <- merge(row_data ,gene_symb, by = "ENTREZID", all.x=TRUE)

Abbiamo poi costruito i rowRanges cioè le posizioni dei geni sul cromosoma di riferimento e il cromosoma di riferimento, per fare ciò abbiamo sfruttato il pacchetto biomaRt che ci permette di collegarci ad un database di Ensembl in particolare a quello del topo, successivamente abbiamo effettuato il mapping grazie al pacchetto GenomicRanges. Infine è stato possibile costruire il summarized Experiment contente tutte le informazioni neccessarie per l’analisi, andiamo poi a salvare il SummarizedExperiment in modo da poterlo richiamare facilmente nelle prossime occasioni

library(dplyr)
library(biomaRt)
library(GenomicRanges)
fonte <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
posizioni_geni <- getBM(filters = "entrezgene_id", attributes = c("entrezgene_id", "chromosome_name", "start_position", "end_position"), values = id_geni, mart = fonte)
names(posizioni_geni)[1] <- "ENTREZID"

#nel caso di duplicati teniamo in considerazione solo la prima occorrenza
posizioni_geni <- posizioni_geni %>% distinct(ENTREZID, .keep_all = TRUE)
row_data <- merge(row_data, posizioni_geni, by = "ENTREZID", all.x=TRUE)
row_data <- row_data[order(row_data$ordine),]
row_data$ordine <- NULL
silenziamento <- lapply(new_names, function(name){
  if (substr(name, 1, 4) == "BRCA") {
    return("S")
  } else {
    return("NS")
  }
})

silenziamento <- unlist(silenziamento)
colData <- cbind(new_names, silenziamento)

se <- SummarizedExperiment(
    assays = conteggi,
    rowData = row_data,
    colData = colData,
)
head(assay(se))
##           BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2 WT-3 WT-4 WT-5
## 100009600      23      29      57      28      28   16   24   23   14   10
## 100009614       0       0       0       0       0    0    0    0    0    0
## 100009664       1       1       0       0       0    0    0    0    0    0
## 100017       2330    2586    2683    2705    3069 1659 1978 1973 1437 1591
## 100019       3428    3522    2736    3616    3941 2353 3357 2990 2366 3362
## 100033459     217     265     166     101     119   36  153   77   59   52
colData(se)
## DataFrame with 10 rows and 2 columns
##           new_names silenziamento
##         <character>   <character>
## BRCA2-1     BRCA2-1             S
## BRCA2-2     BRCA2-2             S
## BRCA2-3     BRCA2-3             S
## BRCA2-4     BRCA2-4             S
## BRCA2-5     BRCA2-5             S
## WT-1           WT-1            NS
## WT-2           WT-2            NS
## WT-3           WT-3            NS
## WT-4           WT-4            NS
## WT-5           WT-5            NS
saveRDS(se, file = "GSE137818_SummarizedExperiment.rds")
rowData(se)
## DataFrame with 20637 rows and 5 columns
##            ENTREZID        SYMBOL chromosome_name start_position end_position
##           <integer>   <character>     <character>      <integer>    <integer>
## 100009600 100009600         Zglp1               9       20973689     20978389
## 100009614 100009614       Gm10024              10       77547291     77547843
## 100009664 100009664 F630206G17Rik              NA             NA           NA
## 100017       100017       Ldlrap1               4      134468865    134495335
## 100019       100019          Mdn1               4       32657119     32775217
## ...             ...           ...             ...            ...          ...
## 99889         99889        Arfip1               3       84403400     84489932
## 99890         99890         Prmt6               3      110153425    110158314
## 99899         99899         Ifi44               3      151436559    151455597
## 99929         99929        Tiparp               3       65435831     65462939
## 99982         99982         Kdm1a               4      136277851    136330034
sort(table(rowData(se)$chromosome_name))
## 
## GL456221.1          Y         18         16         19         12          X 
##          1          2        521        636        660        683        703 
##         14         13         15         10         17          8          3 
##        711        777        805        954        988       1008       1010 
##          6          9          1          5          4          2          7 
##       1050       1101       1190       1208       1234       1556       1580 
##         11 
##       1612
se <- readRDS("GSE137818_SummarizedExperiment.rds")

Step 3: Analisi esplorativa dei dati e filtraggio

Possiamo fare un primo boxplot su scala logaritmica per capire se i dati che guardiamo sono comparabili

library(RColorBrewer) 
pal<-brewer.pal(3, "Set2")
boxplot(log1p(assay(se)), col=pal[as.factor(se$silenziamento)], las=2)

e osserviamo che tutto sommato per tutti i campioni siamo sugli stessi ordini di grandezza, possiamo di seguito effettuare un filtraggio andando a scartare i geni che sono poco presenti in particolare quelli con somma di riga minore o uguale a 10

sef <- se[rowSums(assay(se))>10,]
assay(sef,"log") <- log1p(assay(sef)) 
boxplot(assay(sef,"log"), col=pal[as.factor(sef$silenziamento)], las=2)

A seguito del filtraggio possiamo osservare che rimuovendo i geni scarsamente espressi che potrebbero rappresentare un rumore nei nostri dati, abbiamo migliorato le distribuzioni dei dati perchè è come se avessimo delle distribuzioni più “regolari” e simmetriche rispetto al boxplot precedente

Come un boxplot va interpretato in termini di distribuzione di probabilità
Come un boxplot va interpretato in termini di distribuzione di probabilità

Andiamo poi ad fare un grafico RLE (Relative Log Expression) che ci restituisce come informazione il logaritmo del valore dei nostri geni normalizzati rispetto alla loro mediana su tutti i campioni (log(gi_c1/gi_cmedian)), possiamo osservare che già senza nessuna normalizzazione dovremmo avere una lieve sovraespressione nei campioni con silenziamento di BRCA2 a meno di differenze legate all’atto dell’esperimento, effetto che proveremo poi a mitigare senza appiattire troppo il contenuto informativo del segnale successivamente con la normalizzazione.

library(EDASeq)
plotRLE(as.matrix(assay(sef)), col=pal[as.factor(sef$silenziamento)], outline=FALSE, las=2)

Andando poi a costruire la PCA osserviamo la netta separazione delle due tipologie di campione, anche se i campioni WT risultano poco vicini in termini della PC2 e la varianza totale spiegata è di circa 50% rispetto al 70-80% che ci aspettiamo per una buona rappresentazione dei dati

plotPCA(as.matrix(assay(sef)), col=pal[as.factor(sef$silenziamento)])

Risulta quindi necessario a valle di queste ultime osservazioni cercare di migliorare la situazione andando a normalizzare i dati

Step 4: Esplorazione delle normalizzazioni

Calcoliamo diverse tipologie di normalizzazione ovvero:

-Upper Quartile: Allinea i dati rispetto al valore del 3 quartile

-Full Quartile: Usa tutti i quartili per calcolare un fattore di scala

-Trimmed Mean of M values: Calcola un fattore di scala sfruttando la media dei rapporti di espressione dei geni

-RLE: “median library is calculated from the geometric mean of all columns and the median ratio of each sample to the median library is taken as the scale factor.”

library(edgeR)
assayNames(sef)[1] <- "counts"
assay(sef, "upper") <- betweenLaneNormalization(as.matrix(assay(sef)), which="upper")
assay(sef, "full") <- betweenLaneNormalization(as.matrix(assay(sef)), which="full")
assay(sef, "tmm") <- cpm(assay(sef), lib.size = calcNormFactors(assay(sef), method = "TMM") * colSums(assay(sef)))
assay(sef, "rle_gm") <- cpm(assay(sef), lib.size = calcNormFactors(assay(sef), method = "RLE") * colSums(assay(sef)))
for(a in c("counts","upper", "full", "tmm", "rle_gm")) {
 plotRLE(as.matrix(assay(sef, a)), col=pal[as.factor(sef$silenziamento)], las=2, outline=FALSE, main = a)
}

for(a in c("counts","upper", "full", "tmm", "rle_gm")) {
 plotPCA(as.matrix(assay(sef, a)), col=pal[as.factor(sef$silenziamento)], main = a)
}

Step5: Stima e validazione dei parametri di normalizzazione e dispersione

A seguito dello step precedente, la normalizzazione che spiega la maggior parte della varianza, permette una netta separazione dei campioni in due regioni dello spazio delle componenti principali senza far collassare i diversi campioni in un unico punto e allo stesso tempo permette nei boxplot di avere mediane allineate e box confrontabili è la “TMM”. Andiamo quindi a calcolare i fattori di normalizzazione e li salviamo nella variabile ‘dge’

dge <- calcNormFactors(sef, method = "TMM")

Andiamo poi a costruire la model matrix, cioè come le covariate (cioè le variabili indipendenti che nel nostro caso è il Silenziamento) sono organizzate (se più di 1, non come nel nostro caso) per modellare la loro dipendenza dalla variabile dipendente che nel nostro caso è l’espressione genica dei vari geni

design <- model.matrix(~ silenziamento ,data=colData(sef))

Andiamo poi a stimare la dispersione dei dati grazie alla funzione “estimateDisp” che sfruttando i fattori di normalizzazione e la matrice di design ci restituisce:

-la dispersione globale comune a tutti i geni

-relazione tra media e dispersione di ogni singolo gene (come la varianza cambia in funzione del livello medio di espressione genica)

-la dispersione specifica per ogni gene

Successivamente andiamo a plottare il grafico “Media-Varianza” in cui:

-Gli elementi grigi rappresentano la varianza grezza cioè la varianza osservata per ciascun gene

-Gli elementi azzurri rappresentano la varianza stimata per ciascun gene sfruttando i parametri di dispersione calcolati nello stesso blocco di codice

-La trend-line rappresenta la varianza attesa dato un certo livello medio di espressione genica, considerando una distribuzione binomiale negativa

dge <- estimateDisp(dge, design)
plotMeanVar(dge, show.raw.vars = TRUE, show.tagwise.vars = TRUE,
 show.ave.raw.vars = FALSE)

Possiamo quindi dedurre dal grafico Media-Varianza che sia i dati grezzi che quelli stimati presentano una varianza maggiore del normale (rispetto alla linea nera), suggerendo la presenza di geni con variazioni significative di espressione tra i vari campioni che potrebbero quindi essere differenzialmente espressi.

A seguire, andiamo a plottare il coefficiente di variazione biologica (BCV) che plotta i conteggi rispetto ad un parametro che misura la dispersione rispetto alla variabilità biologica dei campioni Nel grafico possiamo osservare:

-I pallini neri che rappresentano il coefficiente di variazione biologica stimata per ogni campione

-il coefficiente di variazione biologica comune stimato su tutti i geni

-il coefficiente di variazione biologica stimata rispetto al livello medio di espressione di tutti i geni

plotBCV(dge)

Possiamo a valle di questi grafici fare le seguenti osservazioni: Osserviamo che i tagwise seguono l’andamento della Trend per la maggior parte dei geni, inoltre osserviamo che seppur nel grafico Mean-Variance precedente i geni non fossero perfettamente distribuiti intorno alla varianza attesa, il trend viene seguito in maniera abbastanza fedele, dunque queste informazioni ci danno sicurezza rispetto alla “validazione” dei parametri di stima calcolati con la funzione ‘estimateDisp’

Step 6: Analisi dell’espressione differenziale

Con i parametri di dispersione (per catturare la variabilità biologica) e i coefficienti di normalizzazione (per eliminare variazioni non biologiche) ottenuti e validati nello step precedente, andiamo a costruire lo specifico modello lineare generalizzato (GLM) specificando il contrasto per indicare che nel calcolo del log(FC) debba considerare Silenzati/Non Silenziati e cambiando la matrice di design aggiungendo uno 0 in modo da usare un “non-intercept model” come di consueto in questo tipo di analisi.

design0 <- model.matrix(~ 0 + silenziamento ,data=colData(sef))
cont <- makeContrasts(silenziamentoS - silenziamentoNS, levels=design0)
glm_s <- glmFit(dge, design0)

Eseguiamo il test di espressione differenziale utilizzando il modello costruito

ris <- glmLRT(glm_s, contrast = cont)
top <- topTags(ris, n=Inf)$table
diff_exp <- top[top$FDR<=0.05,]
up_reg <- diff_exp[(diff_exp$logFC)>0,]
down_reg <- diff_exp[(diff_exp$logFC)<0,]
table(top$FDR<=0.05)
## 
## FALSE  TRUE 
## 12872  4273

Ottenendo 4273 geni differenzialmente espressi considerando come p-value ajusted (FDR) ≤ 0.05. che si distribuiscono sui diversi cromosomi nel seguente modo:

sort(table(diff_exp$chromosome_name))
## 
## GL456221.1         18         16         19         12         13          X 
##          1         89        129        137        138        139        141 
##         14         17          9          6          8         10         15 
##        143        185        202        203        204        205        214 
##          4          5          2          3          7         11          1 
##        243        263        289        297        300        327        341

Inoltre abbiamo suddiviso il nostro set in due partizioni, i geni up-regolati che rappresentano una maggiore presenza nel caso di silenziamento di BRCA2 e i geni down-regolati che rappresentano una minore presenza nel caso di silenziamento di BRCA2

dim(up_reg)
## [1] 2276   10
dim(down_reg)
## [1] 1997   10

Step 7: Analisi funzionale ed enrchiment

Ottenuti i geni che risultano espressi in maniera significativamente diversa tra i due cluster di campioni, possiamo andare a fare un’analisi funzionale per cercare di capire a che servono questi geni e quindi in che modo la loro diversa espressione si ripercuote poi in termini funzionali

Analisi Generale

La prima cosa è andare a guardare in generale, senza distinzione di up e down regolati cosa possiamo dedurre

Gene Ontology

Andiamo a fare un’analisi di arricchimento rispetto alla Gene Ontology, e successivamente andiamo a fare diversi plot in diverse condizioni considerando per il cnetplot una partizione relativa ai primi 15 elementi relativi alla gene ontology per avere una maggiore leggibilità

library(simplifyEnrichment)
library(clusterProfiler)
library(enrichplot)

gene_sym <- diff_exp$SYMBOL
go_sym <- enrichGO(gene_sym, OrgDb = organism,
                           keyType ="SYMBOL", ont = "ALL")
dotplot(go_sym, showCategory=10)

go_sym1<-pairwise_termsim(go_sym)
emapplot(go_sym1, showCategory=10)

gene_list <- diff_exp
fold_changes <- setNames(gene_list$logFC, gene_list$SYMBOL)
go_sym_sub <- subset_enrichResult(go_sym, 15)
cnetplot(go_sym, foldChange=fold_changes, circular = TRUE, colorEdge = TRUE)
cnetplot(go_sym_sub, foldChange=fold_changes, circular = TRUE, colorEdge = TRUE)

Usando invece i plot offerti da gprofiler2

library(gprofiler2)
geneOnt_ris <- gost(gene_sym, sources="GO",organism = "mmusculus", significant = TRUE, user_threshold = 0.05)
gostplot(geneOnt_ris, interactive = TRUE)

Andando però a dividere le tre macro aree del vocabolario gene ontology osserviamo che

table(go_sym@result$ONTOLOGY)
## 
##   BP   CC   MF 
## 2607  186  312
((go_sym@result)[go_sym@result$ONTOLOGY == "BP", ])[1:5,3]
## [1] "regulation of innate immune response"              
## [2] "positive regulation of innate immune response"     
## [3] "immune response-regulating signaling pathway"      
## [4] "positive regulation of response to biotic stimulus"
## [5] "response to virus"
((go_sym@result)[go_sym@result$ONTOLOGY == "CC", ])[1:5,3]
## [1] "collagen-containing extracellular matrix"
## [2] "cell-substrate junction"                 
## [3] "focal adhesion"                          
## [4] "cell leading edge"                       
## [5] "membrane raft"
((go_sym@result)[go_sym@result$ONTOLOGY == "MF", ])[1:5,3]
## [1] "GTPase regulator activity"                   
## [2] "nucleoside-triphosphatase regulator activity"
## [3] "GTP binding"                                 
## [4] "guanyl nucleotide binding"                   
## [5] "guanyl ribonucleotide binding"

KEGG Pathways - Overview

Sfruttiamo poi i KEGG Pathways per continuare la nostra analisi fuzionale usando sia ClusterProfiler che Gprofiler2

gene_entrez <- diff_exp$ENTREZID
Enrich.KEG <- enrichKEGG(gene_entrez, organism="mmu")
Enrich.KEG1<-pairwise_termsim(Enrich.KEG)
dotplot(Enrich.KEG , showCategory=10)
emapplot(Enrich.KEG1, , showCategory=10)
cnetplot(Enrich.KEG, foldChange=fold_changes, circular = TRUE, colorEdge = TRUE) 

Kegg_res <- gost(gene_entrez, sources="KEGG", organism = "mmusculus", significant = TRUE,user_threshold = 0.05)
gostplot(Kegg_res, interactive = TRUE)

Abbiamo trovato quindi i seguenti risultati per Gprofiler2 e ClusterProfiler

Kegg_res$result$term_name
##  [1] "Toxoplasmosis"                                         
##  [2] "Focal adhesion"                                        
##  [3] "Th17 cell differentiation"                             
##  [4] "T cell receptor signaling pathway"                     
##  [5] "Th1 and Th2 cell differentiation"                      
##  [6] "Fluid shear stress and atherosclerosis"                
##  [7] "ECM-receptor interaction"                              
##  [8] "MAPK signaling pathway"                                
##  [9] "Fc gamma R-mediated phagocytosis"                      
## [10] "Rheumatoid arthritis"                                  
## [11] "Yersinia infection"                                    
## [12] "Chemokine signaling pathway"                           
## [13] "Leishmaniasis"                                         
## [14] "TNF signaling pathway"                                 
## [15] "Coronavirus disease - COVID-19"                        
## [16] "Pertussis"                                             
## [17] "Efferocytosis"                                         
## [18] "Human papillomavirus infection"                        
## [19] "Measles"                                               
## [20] "Proteoglycans in cancer"                               
## [21] "PI3K-Akt signaling pathway"                            
## [22] "Chagas disease"                                        
## [23] "Influenza A"                                           
## [24] "NOD-like receptor signaling pathway"                   
## [25] "Choline metabolism in cancer"                          
## [26] "PD-L1 expression and PD-1 checkpoint pathway in cancer"
## [27] "Lipid and atherosclerosis"                             
## [28] "Toll-like receptor signaling pathway"                  
## [29] "C-type lectin receptor signaling pathway"              
## [30] "Tuberculosis"                                          
## [31] "Osteoclast differentiation"                            
## [32] "Axon guidance"                                         
## [33] "Metabolic pathways"                                    
## [34] "Cytokine-cytokine receptor interaction"                
## [35] "Leukocyte transendothelial migration"                  
## [36] "Sphingolipid signaling pathway"                        
## [37] "Fc epsilon RI signaling pathway"                       
## [38] "Lysosome"                                              
## [39] "Virion - Herpesvirus"                                  
## [40] "Rap1 signaling pathway"                                
## [41] "Inflammatory bowel disease"                            
## [42] "Arrhythmogenic right ventricular cardiomyopathy"       
## [43] "Amoebiasis"                                            
## [44] "Regulation of actin cytoskeleton"                      
## [45] "Ras signaling pathway"                                 
## [46] "Endocrine resistance"                                  
## [47] "Pathways in cancer"
Enrich.KEG$Description
##  [1] "Toxoplasmosis - Mus musculus (house mouse)"                                                
##  [2] "Focal adhesion - Mus musculus (house mouse)"                                               
##  [3] "Th17 cell differentiation - Mus musculus (house mouse)"                                    
##  [4] "T cell receptor signaling pathway - Mus musculus (house mouse)"                            
##  [5] "Th1 and Th2 cell differentiation - Mus musculus (house mouse)"                             
##  [6] "MAPK signaling pathway - Mus musculus (house mouse)"                                       
##  [7] "Fluid shear stress and atherosclerosis - Mus musculus (house mouse)"                       
##  [8] "ECM-receptor interaction - Mus musculus (house mouse)"                                     
##  [9] "Fc gamma R-mediated phagocytosis - Mus musculus (house mouse)"                             
## [10] "Rheumatoid arthritis - Mus musculus (house mouse)"                                         
## [11] "Chemokine signaling pathway - Mus musculus (house mouse)"                                  
## [12] "TNF signaling pathway - Mus musculus (house mouse)"                                        
## [13] "Yersinia infection - Mus musculus (house mouse)"                                           
## [14] "Leishmaniasis - Mus musculus (house mouse)"                                                
## [15] "Measles - Mus musculus (house mouse)"                                                      
## [16] "Efferocytosis - Mus musculus (house mouse)"                                                
## [17] "Human papillomavirus infection - Mus musculus (house mouse)"                               
## [18] "Pertussis - Mus musculus (house mouse)"                                                    
## [19] "Proteoglycans in cancer - Mus musculus (house mouse)"                                      
## [20] "PI3K-Akt signaling pathway - Mus musculus (house mouse)"                                   
## [21] "Chagas disease - Mus musculus (house mouse)"                                               
## [22] "Influenza A - Mus musculus (house mouse)"                                                  
## [23] "Cytoskeleton in muscle cells - Mus musculus (house mouse)"                                 
## [24] "Toll-like receptor signaling pathway - Mus musculus (house mouse)"                         
## [25] "Choline metabolism in cancer - Mus musculus (house mouse)"                                 
## [26] "PD-L1 expression and PD-1 checkpoint pathway in cancer - Mus musculus (house mouse)"       
## [27] "Lipid and atherosclerosis - Mus musculus (house mouse)"                                    
## [28] "Tuberculosis - Mus musculus (house mouse)"                                                 
## [29] "NOD-like receptor signaling pathway - Mus musculus (house mouse)"                          
## [30] "Cytokine-cytokine receptor interaction - Mus musculus (house mouse)"                       
## [31] "C-type lectin receptor signaling pathway - Mus musculus (house mouse)"                     
## [32] "Axon guidance - Mus musculus (house mouse)"                                                
## [33] "Osteoclast differentiation - Mus musculus (house mouse)"                                   
## [34] "Sphingolipid signaling pathway - Mus musculus (house mouse)"                               
## [35] "Leukocyte transendothelial migration - Mus musculus (house mouse)"                         
## [36] "Lysosome - Mus musculus (house mouse)"                                                     
## [37] "Fc epsilon RI signaling pathway - Mus musculus (house mouse)"                              
## [38] "Virion - Herpesvirus - Mus musculus (house mouse)"                                         
## [39] "Rap1 signaling pathway - Mus musculus (house mouse)"                                       
## [40] "Inflammatory bowel disease - Mus musculus (house mouse)"                                   
## [41] "Ras signaling pathway - Mus musculus (house mouse)"                                        
## [42] "Regulation of actin cytoskeleton - Mus musculus (house mouse)"                             
## [43] "Amoebiasis - Mus musculus (house mouse)"                                                   
## [44] "Endocrine resistance - Mus musculus (house mouse)"                                         
## [45] "Arrhythmogenic right ventricular cardiomyopathy - Mus musculus (house mouse)"              
## [46] "Alcoholic liver disease - Mus musculus (house mouse)"                                      
## [47] "Central carbon metabolism in cancer - Mus musculus (house mouse)"                          
## [48] "B cell receptor signaling pathway - Mus musculus (house mouse)"                            
## [49] "HIF-1 signaling pathway - Mus musculus (house mouse)"                                      
## [50] "Epstein-Barr virus infection - Mus musculus (house mouse)"                                 
## [51] "Neurotrophin signaling pathway - Mus musculus (house mouse)"                               
## [52] "Primary immunodeficiency - Mus musculus (house mouse)"                                     
## [53] "Glycerophospholipid metabolism - Mus musculus (house mouse)"                               
## [54] "RIG-I-like receptor signaling pathway - Mus musculus (house mouse)"                        
## [55] "Viral protein interaction with cytokine and cytokine receptor - Mus musculus (house mouse)"
## [56] "ErbB signaling pathway - Mus musculus (house mouse)"                                       
## [57] "mTOR signaling pathway - Mus musculus (house mouse)"                                       
## [58] "Asthma - Mus musculus (house mouse)"                                                       
## [59] "Cytosolic DNA-sensing pathway - Mus musculus (house mouse)"                                
## [60] "Protein processing in endoplasmic reticulum - Mus musculus (house mouse)"                  
## [61] "Estrogen signaling pathway - Mus musculus (house mouse)"                                   
## [62] "Hematopoietic cell lineage - Mus musculus (house mouse)"                                   
## [63] "Relaxin signaling pathway - Mus musculus (house mouse)"                                    
## [64] "Apoptosis - Mus musculus (house mouse)"                                                    
## [65] "Malaria - Mus musculus (house mouse)"                                                      
## [66] "Adipocytokine signaling pathway - Mus musculus (house mouse)"                              
## [67] "Thyroid hormone signaling pathway - Mus musculus (house mouse)"                            
## [68] "Prostate cancer - Mus musculus (house mouse)"                                              
## [69] "Hypertrophic cardiomyopathy - Mus musculus (house mouse)"                                  
## [70] "Cell adhesion molecules - Mus musculus (house mouse)"                                      
## [71] "AGE-RAGE signaling pathway in diabetic complications - Mus musculus (house mouse)"         
## [72] "Hepatitis C - Mus musculus (house mouse)"                                                  
## [73] "Pyruvate metabolism - Mus musculus (house mouse)"                                          
## [74] "Sphingolipid metabolism - Mus musculus (house mouse)"                                      
## [75] "Salmonella infection - Mus musculus (house mouse)"                                         
## [76] "JAK-STAT signaling pathway - Mus musculus (house mouse)"                                   
## [77] "IL-17 signaling pathway - Mus musculus (house mouse)"                                      
## [78] "Hepatitis B - Mus musculus (house mouse)"                                                  
## [79] "Tight junction - Mus musculus (house mouse)"                                               
## [80] "VEGF signaling pathway - Mus musculus (house mouse)"                                       
## [81] "Endocytosis - Mus musculus (house mouse)"                                                  
## [82] "NF-kappa B signaling pathway - Mus musculus (house mouse)"                                 
## [83] "Phospholipase D signaling pathway - Mus musculus (house mouse)"                            
## [84] "Dilated cardiomyopathy - Mus musculus (house mouse)"                                       
## [85] "Legionellosis - Mus musculus (house mouse)"                                                
## [86] "Platinum drug resistance - Mus musculus (house mouse)"                                     
## [87] "Insulin resistance - Mus musculus (house mouse)"                                           
## [88] "Kaposi sarcoma-associated herpesvirus infection - Mus musculus (house mouse)"              
## [89] "Synaptic vesicle cycle - Mus musculus (house mouse)"                                       
## [90] "Glycosaminoglycan biosynthesis - keratan sulfate - Mus musculus (house mouse)"             
## [91] "Hepatocellular carcinoma - Mus musculus (house mouse)"                                     
## [92] "Small cell lung cancer - Mus musculus (house mouse)"                                       
## [93] "Coronavirus disease - COVID-19 - Mus musculus (house mouse)"                               
## [94] "African trypanosomiasis - Mus musculus (house mouse)"                                      
## [95] "Antigen processing and presentation - Mus musculus (house mouse)"                          
## [96] "Biosynthesis of amino acids - Mus musculus (house mouse)"                                  
## [97] "Human immunodeficiency virus 1 infection - Mus musculus (house mouse)"

Che sono coerenti tra loro, possiamo inoltre notare che come visto a lezione è presente il pathway per il Covid19 ma non possiamo dire se questo silenziamento in termini di KEGG Pathways vada ad influire positivamente o nega tivamente per il cancro al seno, infatti come possiamo notare è presente “Pathways in cancer” ([47] di Kegg_res) ma non sappiamo se in termini positivi o negativi

KEGG Pathways - Geni UP REGOLATI

upgene_entrez <- up_reg$ENTREZID
Kegg_res <- gost(upgene_entrez, sources="KEGG", organism = "mmusculus", significant = TRUE,user_threshold = 0.05)
Kegg_res$result$term_name
##  [1] "Th1 and Th2 cell differentiation"                             
##  [2] "Th17 cell differentiation"                                    
##  [3] "Toxoplasmosis"                                                
##  [4] "NOD-like receptor signaling pathway"                          
##  [5] "Influenza A"                                                  
##  [6] "Measles"                                                      
##  [7] "Epstein-Barr virus infection"                                 
##  [8] "T cell receptor signaling pathway"                            
##  [9] "Leishmaniasis"                                                
## [10] "Coronavirus disease - COVID-19"                               
## [11] "Inflammatory bowel disease"                                   
## [12] "Staphylococcus aureus infection"                              
## [13] "Herpes simplex virus 1 infection"                             
## [14] "Cell adhesion molecules"                                      
## [15] "Pertussis"                                                    
## [16] "Cytokine-cytokine receptor interaction"                       
## [17] "Chemokine signaling pathway"                                  
## [18] "Cytosolic DNA-sensing pathway"                                
## [19] "Toll-like receptor signaling pathway"                         
## [20] "Rheumatoid arthritis"                                         
## [21] "Osteoclast differentiation"                                   
## [22] "Asthma"                                                       
## [23] "Yersinia infection"                                           
## [24] "Chagas disease"                                               
## [25] "Primary immunodeficiency"                                     
## [26] "Fc gamma R-mediated phagocytosis"                             
## [27] "Tuberculosis"                                                 
## [28] "Intestinal immune network for IgA production"                 
## [29] "RIG-I-like receptor signaling pathway"                        
## [30] "Hepatitis C"                                                  
## [31] "Leukocyte transendothelial migration"                         
## [32] "TNF signaling pathway"                                        
## [33] "JAK-STAT signaling pathway"                                   
## [34] "Virion - Herpesvirus"                                         
## [35] "Fc epsilon RI signaling pathway"                              
## [36] "Viral protein interaction with cytokine and cytokine receptor"
## [37] "PD-L1 expression and PD-1 checkpoint pathway in cancer"       
## [38] "Antigen processing and presentation"

KEGG Pathways - Geni DOWN REGOLATI

downgene_entrez <- down_reg$ENTREZID
Kegg_res <- gost(downgene_entrez, sources="KEGG", organism = "mmusculus", significant = TRUE,user_threshold = 0.05)
Kegg_res$result$term_name
##  [1] "Focal adhesion"                             
##  [2] "HIF-1 signaling pathway"                    
##  [3] "MAPK signaling pathway"                     
##  [4] "Protein processing in endoplasmic reticulum"
##  [5] "Fluid shear stress and atherosclerosis"     
##  [6] "PI3K-Akt signaling pathway"                 
##  [7] "mTOR signaling pathway"                     
##  [8] "Proteoglycans in cancer"                    
##  [9] "Lysosome"                                   
## [10] "Central carbon metabolism in cancer"        
## [11] "Metabolic pathways"                         
## [12] "Choline metabolism in cancer"               
## [13] "ECM-receptor interaction"                   
## [14] "Sphingolipid signaling pathway"             
## [15] "Rap1 signaling pathway"                     
## [16] "Autophagy - animal"                         
## [17] "Efferocytosis"                              
## [18] "Thyroid hormone signaling pathway"          
## [19] "Renal cell carcinoma"                       
## [20] "Amoebiasis"

KEGG Pathways - Considerazioni finali

“La sovraespressione di pathways immunitari e infiammatori suggerisce una maggiore attivazione della risposta immunitaria, che può essere positiva se il sistema immunitario riesce a riconoscere e attaccare le cellule tumorali. Tuttavia, l’infiammazione cronica può anche favorire la progressione tumorale.

La sottoespressione di pathways legati alla crescita, adesione e metabolismo cellulare suggerisce una ridotta capacità del tumore di proliferare, sopravvivere e metastatizzare, il che è generalmente positivo per il controllo del tumore. Conclusione:

Nel complesso, il silenziamento di BRCA2 sembra avere un effetto misto, ma prevalentemente positivo rispetto al controllo del tumore al seno. La riduzione della proliferazione e della capacità metastatica delle cellule tumorali, insieme a una potenziale maggiore suscettibilità alle terapie immunitarie, suggerisce che il silenziamento di BRCA2 potrebbe limitare la progressione del tumore al seno. Tuttavia, è essenziale considerare il contesto specifico e ulteriori studi sperimentali per confermare questa valutazione.”

Un’ulteriore informazione interessante è che separando up e down regolati otteniamo un numero minore di pathway rispetto al considerarli entrambi contemporaneamente

#install.packages("remotes")
#remotes::install_github("jokergoo/simplifyEnrichment")
#per poter usare la funzione che fa un subset di un enrichResult: subset_enrichResult(elemento, quantità)

Step 8: Network Analysis

Similarità dei campioni

Come ultima parte di questo progetto, ci si pone l’obbiettivo di analizzare i nostri campioni in termini di rete per vedere come questi di distribuiscono in termini di similarità, prima di tutto andiamo calcolare la matrice delle distanze della trasposta dei nostri geni

library(igraph)
library(philentropy)
Gene_Expression <- assay(sef, "tmm")
dist <- as.matrix(distance(t(Gene_Expression), method = "euclidean", use.row.names = TRUE))

ottendo le seguenti distanze

dist
##           BRCA2-1   BRCA2-2   BRCA2-3   BRCA2-4   BRCA2-5      WT-1      WT-2
## BRCA2-1     0.000  3950.972  4078.287  5861.096  3395.781 16882.499 11319.844
## BRCA2-2  3950.972     0.000  3082.494  7578.420  4259.664 20038.013 13630.750
## BRCA2-3  4078.287  3082.494     0.000  7603.953  4353.241 19291.317 13088.773
## BRCA2-4  5861.096  7578.420  7603.953     0.000  4801.625 15145.999  9625.304
## BRCA2-5  3395.781  4259.664  4353.241  4801.625     0.000 17161.324 10957.973
## WT-1    16882.499 20038.013 19291.317 15145.999 17161.324     0.000  9676.455
## WT-2    11319.844 13630.750 13088.773  9625.304 10957.973  9676.455     0.000
## WT-3    15141.359 17802.041 17175.606 13476.008 15245.875  6690.066  5241.374
## WT-4    20013.154 22873.268 22292.685 17756.943 20085.051  6072.720 10596.133
## WT-5    20258.509 23139.703 22624.269 17262.267 20062.953  6839.825 11533.335
##              WT-3      WT-4      WT-5
## BRCA2-1 15141.359 20013.154 20258.509
## BRCA2-2 17802.041 22873.268 23139.703
## BRCA2-3 17175.606 22292.685 22624.269
## BRCA2-4 13476.008 17756.943 17262.267
## BRCA2-5 15245.875 20085.051 20062.953
## WT-1     6690.066  6072.720  6839.825
## WT-2     5241.374 10596.133 11533.335
## WT-3        0.000  6295.002  8437.072
## WT-4     6295.002     0.000  4684.591
## WT-5     8437.072  4684.591     0.000

Rispetto alle distanze ottenute, andiamo a “tagliare” i collegamenti per quelli che cadono al di sotto del quantile calcolato al 50% per non avere una rete densa e per poter identificare i sottogruppi, inoltre la nostra rete non è pesata quindi sostituiamo al valore di distanza 1 per indicare la presenza di un collegamento ed inoltre annulliamo anche la diagonale principale per evitare loop locali

dist[dist<quantile(dist,0.5)] <- 1
dist[dist>=quantile(dist,0.5)] <- 0
diag(dist) <- 0
dist
##         BRCA2-1 BRCA2-2 BRCA2-3 BRCA2-4 BRCA2-5 WT-1 WT-2 WT-3 WT-4 WT-5
## BRCA2-1       0       1       1       1       1    0    0    0    0    0
## BRCA2-2       1       0       1       1       1    0    0    0    0    0
## BRCA2-3       1       1       0       1       1    0    0    0    0    0
## BRCA2-4       1       1       1       0       1    0    1    0    0    0
## BRCA2-5       1       1       1       1       0    0    0    0    0    0
## WT-1          0       0       0       0       0    0    1    1    1    1
## WT-2          0       0       0       1       0    1    0    1    1    0
## WT-3          0       0       0       0       0    1    1    0    1    1
## WT-4          0       0       0       0       0    1    1    1    0    1
## WT-5          0       0       0       0       0    1    0    1    1    0

a questo punto abbiamo una matrice che ci permette di costruire un grafico in maniera leggera in termini computazionali costruendo un oggeto di tipo iGraph

net <- graph_from_adjacency_matrix(dist, mode="undirected")
plot(net, vertex.size=5,vertex.label.cex=0.5, vertex.frame.color="#ffffff",vertex.color="green")

Osservando come si creino due reti ben separate tra i campioni silenziati e quelli non silenziati dando validità in termini di robustezza alle nostre analisi perchè i campioni dei due gruppi sono simili tra loro.

Clustering dell’enrichmentResult del GO - Ovewview

#remotes::install_github("jmw86069/jamenrich")
library(multienrichjam)
library(jamba)
library(colorjam)
go_net <- (go_sym)[go_sym@result$ONTOLOGY == "BP", ]
go_net<-enrichMapJam(go_net)
plot(go_net, vertex.size=8,vertex.label="", vertex.frame.color="#ffffff",vertex.color="green")

Possiamo osservare come nella nostra rete siano presenti dei gruppi più o meno grandi e vogliamo indentificarli con l’algoritmo Louvaine, previa una rimozione dei nodi isolati

ind <- which((degree(go_net)) == 0)
go_net<- delete_vertices(go_net, ind)
cl <- cluster_louvain(go_net)
plot(go_net, vertex.size=5,vertex.label="", vertex.frame.color="#ffffff",vertex.color=cl$membership) 

Abbiamo quindi identificato i nostri clusters, per capire questi clusters se sono up o down regolati e quindi, nel caso in esame del silenziamento di BRCA2 vanno ad inibire o stimolare determinati processi bisogna considerare l’up e down regolation dei geni e non la semplice espressione differenziale come abbiamo appena fatto

Clustering dell’enrichment Result del GO - Geni UP REGOLATI

upgene_sym <- up_reg$SYMBOL
upgo_sym <- enrichGO(upgene_sym, OrgDb = organism,
                           keyType ="SYMBOL", ont = "ALL")
upgo_net <- (upgo_sym)[upgo_sym@result$ONTOLOGY == "BP", ]
upgo_net<-enrichMapJam(upgo_net)

ind <- which((degree(upgo_net)) == 0)
upgo_net<- delete_vertices(upgo_net, ind)
cl <- cluster_louvain(upgo_net)
plot(upgo_net, vertex.size=5,vertex.label="", vertex.frame.color="#ffffff",vertex.color=cl$membership) 

Identificando 3 cluster ben distinti che possiamo studiare

df <- data.frame(length, index)
for (i in seq_along(cl)) {
  df <- rbind(df, data.frame(length = length(unlist(cl[i])), index = i))
}
df_clust <- df[order(df$length, decreasing = TRUE), ]
df_clust
##   length index
## 2     25     2
## 1     22     1
## 3      2     3
for (i in 1:nrow(df_clust)) {
  num <- df_clust[i, 2]
  cat("Cluster", num, ": è formato da\n")
  cat(paste(cl[[num]], collapse = ",\n"), "\n\n")
}
## Cluster 2 : è formato da
## leukocyte cell-cell adhesion,
## regulation of leukocyte cell-cell adhesion,
## regulation of T cell activation,
## positive regulation of leukocyte cell-cell adhesion,
## lymphocyte differentiation,
## positive regulation of cell-cell adhesion,
## regulation of immune effector process,
## lymphocyte proliferation,
## positive regulation of cell activation,
## leukocyte proliferation,
## mononuclear cell proliferation,
## positive regulation of T cell activation,
## positive regulation of leukocyte activation,
## T cell differentiation,
## T cell proliferation,
## positive regulation of lymphocyte activation,
## regulation of leukocyte differentiation,
## regulation of hemopoiesis,
## regulation of leukocyte proliferation,
## regulation of lymphocyte proliferation,
## regulation of lymphocyte differentiation,
## regulation of mononuclear cell proliferation,
## leukocyte activation involved in immune response,
## positive regulation of immune effector process,
## cell activation involved in immune response 
## 
## Cluster 1 : è formato da
## defense response to virus,
## defense response to symbiont,
## regulation of innate immune response,
## immune response-regulating signaling pathway,
## immune response-activating signaling pathway,
## positive regulation of innate immune response,
## positive regulation of response to biotic stimulus,
## activation of innate immune response,
## response to interferon-beta,
## cellular response to interferon-beta,
## cytokine-mediated signaling pathway,
## innate immune response-activating signaling pathway,
## pattern recognition receptor signaling pathway,
## cytosolic pattern recognition receptor signaling pathway,
## interferon-mediated signaling pathway,
## negative regulation of viral process,
## negative regulation of cytokine production,
## type I interferon production,
## regulation of viral life cycle,
## regulation of inflammatory response,
## antiviral innate immune response,
## response to virus 
## 
## Cluster 3 : è formato da
## peptidyl-tyrosine phosphorylation,
## peptidyl-tyrosine modification

Abbiamo quindi trovato per ora questi 3 cluster di processi biologici che nel caso del silenziamento di BRCA2 sono sovra-stimolati

Clustering dell’enrichment Result del GO- Geni DOWN REGOLATI

downgene_sym <- down_reg$SYMBOL
downgo_sym <- enrichGO(downgene_sym, OrgDb = organism,
                           keyType ="SYMBOL", ont = "ALL")
downgo_net <- (downgo_sym)[downgo_sym@result$ONTOLOGY == "BP", ]
downgo_net<-enrichMapJam(downgo_net)

ind <- which((degree(downgo_net)) == 0)
downgo_net<- delete_vertices(downgo_net, ind)
cl <- cluster_louvain(downgo_net)
plot(downgo_net, vertex.size=5,vertex.label="", vertex.frame.color="#ffffff",vertex.color=cl$membership) 

Identificando 9 cluster ben distinti che possiamo studiare

df <- data.frame(length, index)
for (i in seq_along(cl)) {
  df <- rbind(df, data.frame(length = length(unlist(cl[i])), index = i))
}
df_clust <- df[order(df$length, decreasing = TRUE), ]
df_clust
##   length index
## 1     11     1
## 5      6     5
## 3      4     3
## 4      4     4
## 6      4     6
## 2      3     2
## 7      3     7
## 8      3     8
## 9      3     9
for (i in 1:nrow(df_clust)) {
  num <- df_clust[i, 2]
  cat("Cluster", num, ": è formato da\n")
  cat(paste(cl[[num]], collapse = ",\n"), "\n\n")
}
## Cluster 1 : è formato da
## regulation of vasculature development,
## regulation of angiogenesis,
## epithelial cell migration,
## epithelium migration,
## tissue migration,
## mesenchymal cell differentiation,
## regulation of epithelial cell migration,
## mesenchyme development,
## positive regulation of angiogenesis,
## positive regulation of vasculature development,
## ameboidal-type cell migration 
## 
## Cluster 5 : è formato da
## regulation of actin filament-based process,
## regulation of supramolecular fiber organization,
## actin filament organization,
## regulation of cellular component size,
## regulation of actin cytoskeleton organization,
## regulation of actin filament organization 
## 
## Cluster 3 : è formato da
## ossification,
## bone mineralization,
## biomineral tissue development,
## bone development 
## 
## Cluster 4 : è formato da
## chemotaxis,
## taxis,
## leukocyte migration,
## cell chemotaxis 
## 
## Cluster 6 : è formato da
## response to peptide,
## cellular response to peptide,
## cellular response to peptide hormone stimulus,
## response to peptide hormone 
## 
## Cluster 2 : è formato da
## cell-substrate adhesion,
## cell-matrix adhesion,
## regulation of cell-substrate adhesion 
## 
## Cluster 7 : è formato da
## extracellular matrix organization,
## external encapsulating structure organization,
## extracellular structure organization 
## 
## Cluster 8 : è formato da
## response to endoplasmic reticulum stress,
## response to topologically incorrect protein,
## response to unfolded protein 
## 
## Cluster 9 : è formato da
## positive regulation of cell projection organization,
## regulation of neurogenesis,
## gliogenesis

Abbiamo quindi trovato per ora questi 9 cluster di processi biologici che nel caso del silenziamento di BRCA2 sono inibiti.

Clustering dell’enrichment Result del GO - Considerazioni finali

Alla luce dei risultati possiamo trarre le seguenti conslusioni riguardo all’analisi dell’arricchimento fatta con la gene ontology:

“Nel complesso, i geni up-regolati mostrano un miglioramento della risposta immunitaria che può essere positivo per la sorveglianza e l’eliminazione delle cellule tumorali, ma possono anche promuovere l’infiammazione cronica che può avere effetti pro-tumorali. D’altra parte, i geni down-regolati riducono processi cruciali per la crescita, la sopravvivenza e la metastasi delle cellule tumorali, suggerendo un effetto complessivamente positivo contro il tumore.”

Inoltre:

“La clusterizzazione non è stata inutile; al contrario, ha fornito un quadro più chiaro e dettagliato dei processi biologici coinvolti. Ha permesso di valutare meglio l’impatto del silenziamento di BRCA2 e di identificare pathways specifici che possono essere target terapeutici.

In sintesi, mentre una lista di geni regolati è utile, la clusterizzazione offre un valore aggiunto significativo per l’interpretazione biologica e la valutazione finale.”

Step : Heatmap Analisys

In questo ultimo step, diamo uno sguardo alle heatmap e osserviamo diverse informazioni interessanti

Analisi generale

Guardiamo in primis un po’ di heatmap dei geni differenzialmente espressi in generale, quelli up regolati e down regolati

library(pheatmap)
pheatmap(assay(sef)[rownames(assay(sef,"tmm")) %in% diff_exp$ENTREZID, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, show_colnames = TRUE, 
         main = "Heatmap dei Geni Differenzialmente Espressi")

pheatmap(assay(sef)[rownames(assay(sef,"tmm")) %in% up_reg$ENTREZID, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, show_colnames = TRUE, 
         main = "Heatmap dei Geni Up_regolati")

pheatmap(assay(sef)[rownames(assay(sef,"tmm")) %in% down_reg$ENTREZID, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, show_colnames = TRUE, 
         main = "Heatmap dei Geni Down_regolati")

Osserviamo come ci sia una netta separazione in termini di gruppi tra BRCA e WT ad eccezione della mappa dei down-regolati in cui il campione WT-4 viene clusterizzato nel gruppo dei BRCA mentre il campione BRCA2-5 viene clusterizzato nel gruppo dei WT

Heatmap in un pathway

Dato che si parla di cancro, uno dei pathway sicuramente interessante è “immune response-regulating signaling pathway” identificato nei GO-Terms che risulta up-regolato

upgene_sym <- up_reg$ENTREZID
upgo_sym <- enrichGO(upgene_sym, OrgDb = organism,
                           keyType ="ENTREZID", ont = "BP")
pos<-which(upgo_sym@result$Description == "immune response-regulating signaling pathway")
pathway_genes <- strsplit(upgo_sym@result[pos,]$geneID, split = "/")[[1]]

pheatmap(assay(sef, "tmm")[rownames(assay(sef)) %in% pathway_genes, ], scale = "row", cluster_rows = TRUE, cluster_cols = TRUE,
         show_rownames = FALSE, show_colnames = TRUE, 
         main = "Heatmap dei Geni Up Regolati del Pathway")

Possiamo osservare come, in questo pathway, tutti i geni contribuiscano allo stesso modo perchè non ci sono regioni orizzontali particolamente intense rispetto alle altre (nella partizione BRCA2) anzi, le regioni di maggior intensità per ogni campione di distribuiscono su sotto-gruppi di geni diversi, questo suggerisce che il Pathway di interesse si verifichi effettivamente piuttosto che, come nel caso del pathway del covid 19 questo venga identificato in maniera erronea